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Abstract 



We present a non-perturbative approach to the problem of quasiparticles 
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coupled to spin-fluctuations. If the fully dressed spin-fluctuation propagator is 



used in the Feynman graph expansion of the single particle Green's function, 
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CN ' the problem of summing all spin-fluctuation exchange Feynman graphs (i.e 
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Q\ ' Gaussian distributed random vector fields. The quasiparticle propagator is 
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all the graphs without fermion loops) can be cast as a functional integral over 



then obtained as the average over all field configurations of the non-interacting 
electron Green's function in an external field. A Monte Carlo sampling of 
this functional integral does not suffer from the "fermion sign problem" and 

o" 

offers an attractive alternative to perturbative calculations. We compare the 

results of our computer simulations to perturbation theory and self-consistent 

5— i 

C3 ' Eliashberg calculations 



I. INTRODUCTION 

In previous papers (Monthoux and Pines, 1993, 1994) we described the results of exten- 
sive one-loop self-consistent calculations on the normal state and superconducting transition 
temperature for the optimally doped copper oxide superconductor YBa 2 Cu 3 07. In these cal- 
culations, the effective low-energy quasiparticle-quasiparticle interaction is expressed as the 
exchange of dressed magnetic excitations whose spectrum is determined by fits to NMR 
experiments. These Eliashberg calculations demonstrated that while strong coupling effects 
reduce T c substantially, d x 2_ y 2 superconductivity at high temperatures is still possible with 
a quasiparticle-spin-fluctuation coupling constant g 2 in the range required to explain the 
observed normal-state resisitivity and optical conductivity. 

Questions regarding the validity of the single spin-fluctuation exchange picture have 
been raised (Schrieffer, 1995). The role of vertex corrections in the paramagnetic phase 
has been examined in detail by Chubukov et al. (1997) and the author (1997). For a 
large Fermi surface consistent with ARPES experiments on YBa2Cu307, vertex corrections 
enhance the effective spin-fluctuation interaction as the antiferromagnetic correlation length 
increases. The effect is strongly anisotropic and most important for quasiparticles at the 
hot spots (Hlubina and Rice, 1997) for which the Midgal argument does not carry through. 
A quasiparticle near the hot spot can emit or absorb any number of spin-fluctuations at 
the antiferromagnetic wavevector Q = (71", tc) where the interaction is strong while staying 
close to the Fermi surface. One would thus expect that as the coupling gets stronger or as 
the correlation length increases, two, three, etc.. spin-fluctuation exchange graphs become 
important and perturbation theory breaks down. Non-perturbative methods are thus needed 
to study the behaviour of hot quasiparticles for long correlation lengths. In the case where 
the charactersitic spin-fluctuation frequency lusf is much smaller than 27rT, only the zero 
Matsubara frequency of the spin-fluctuation propagator is relevant and a solution to the 
problem in this quasi-static limit has been proposed by Schmalian et al (1998). 

In this paper, the calculation of the single particle Green's function for an arbitrary 



number of spin-fluctuation exchanges is formulated as the average of the non-interacting 
electron Green's function in an external magnetic field whose probability distribution is 
Gaussian. Upon expanding the Green's function in powers of the magnetic field and doing 
the averages with Wick's theorem, one recovers the standard time-ordered perturbation 
series for the quasiparticle propagator. However, sampling the functional integral by the 
Monte Carlo method combined with a numerical calculation of the Green's function for each 
external field configuration provides an attractive and non-perturbative method of solution 
when straightforward perturbation theory is expected to break down. 

The low-energy effective action for the planar quasiparticles is given by 
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The spin density s(q, r) is given by 

s(q,r) = Yl ^+q,a(^)^, 7 ^p, 7 ( r ) ( 2 ) 
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where a denotes the three Pauli matrices. The quasiparticle dispersion relation is 

e^ = — 2t(cos(p x ) + cos(p y )) - 4t / cos(p x )cos(p y ) (3) 

with hopping matrix elements t=0.25 eV and t'=-0.45t, as suggested by Si et al. (1993) from 
fits to band theory and ARPES experiments. We use units in which the lattice constant is 
unity, fi denotes the chemical potential, j3 the inverse temperature, g 2 the coupling constant 
and ipi and ^ j0 - are Grassmann variables. A coupling constant g 2 = 0.41 eV 2 gives an 
Eliashberg T c to a d x 2_ y 2 superconducting state of 90°K (Monthoux and Pines (1994)). Since 
the interaction x(q, t) is determined from fits to normal state experiments, it is fully dressed. 
Therefore, in a Feynman graph expansion, the diagrams corresponding to a renormalization 
of the spin-fluctuation propagator x(q 3 T ) should not be included. On the real frequency 
axis, %(q, uS) is assumed to take the phenomenological form introduced by Millis et al.(1990) 



X(q, w) = ^ (4) 

l + e(q-Q) 2 -i^ 

Here \Q is the static spin susceptibility at wavevector Q = (7r, it) and £ is a temperature 
dependent antiferromagnetic correlation length. The frequency ujsf defines the charactersitic 
energy of the low- frequency magnetic response. The quantities xq an d cjsf may be expressed 
in terms of the experimentally measured long- wavelength spin susceptibility Xo an d a short- 
distance characteristic spin-fluctuation energy I\f by introducing a scale factor (3 (not to 
be confused with the inverse temperature) 

Xq = XO^C 2 (5) 

USF = 7M> (6) 

We adopt the same values of these parameters as in Monthoux and Pines (1994), namely 
Taf = 1-3 eV, Xo = 2.6 states/eV, j3 = 32. The temperature dependence of the character- 
sitic spin-fluctuation energy ujsf is taken as 

^sf(T) = (9.5 + 4.75[T(°K)/100]) meV (7) 

from which the temperature-dependent antiferromagnetic correlation length can be obtained 
using Eq. (|6|). The numerical values of these parameters for YBa 2 Cu 3 07 are Xq(T c ) ~ 
75 states/eV, f (T c ) w 2.3 and cj S f(T c ) w 14 meV. 

The spin-fluctuation propagator on the imaginary axis, x(q, i^n) is related to the imagi- 
nary part of the response function Imx(q, a;), Eq. (||), via the spectral representation 



/^ . n /+°° dcj Imx(q, cj) 
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To get %(q, iz/ n ) to decay as 1/^ as i/ n — ► oc, as it should, we introduce a cutoff ujq and take 
Im%(q, a;) = for uj > uq. In the following, we adopt the value Uq = 0.4 eV. 



II. FUNCTIONAL INTEGRAL FORMULATION 

In this section, we formulate the problem of calculating the single particle Green's func- 
tion for an arbitrary number of spin-fluctuation exchanges as a functional integral that is 
suitable for computer simulation. 

Let us first introduce some notation. In the following, we denote by x the position and 
imaginary time coordinate, x = (x, r) and similarly, let q stand for the wavevector and 
Matsubara frequency 3- vector q = (q, w n ). Consider a real vector field M(x), such that 

M(-q) = M*(q) (9) 

and introduce the average < ... > over the Gaussian probability distribution of the fields 
M(q): 

<->4//™--(-?^F) <"> 

where a(q) is real and the normalization factor Z is defined as 

The integrals DM are carried out over all of the independent components of the vector field 
M(q). The two-point correlation function of the field M(q) is straightforward to calculate 
and is given by 

. r / x. r /-, x c- c- f 2a(q) if M(q) is complex , . 

< M i ( q )M i (k) >= s { A - Q x \ N .,.;; x . , (12) 

KHJ n ; ' J ' q \ a(q) if M(q) is real v ; 

Since the integrals over M(q) are Gaussian, the average of an odd number of fields vanishes 
and the higher order correlations can be obtained from the two-point correlation function 
by Wick's theorem. For instance the four-point correlation function is given by: 

< M i (q)M j (k)M n (p)M m (l) > = < Mi(q)Mj(k) X M n (p)M m (l) > 

+ < M i (q)M n (p) X Mj(k)M m (l) > 

+ < M i (q)M m (l) X Mj(k)M n (p) > (13) 



Now we choose a(q) such that 

< M^M^-q) >= y^x(q) (14) 

where N is the number of momenta in the Brillouin Zone, T the temperature and let 

G=<G [1-MG ]" 1 > (15) 

where the matrices Go and M in 3-momentum and spin space are defined as 

[G ](q,a)(q',/3) = G (q)5( q , Q ),(q',/3) (16) 

[M] (q , a)(q ^ ) = a ai/3 -M(q-q / ) (17) 

where a are the three Pauli matrices and a, j3 spin indices. Upon expanding and doing the 
averages over the random fields M, one generates the following perturbation expansion: 

G = <Go[l-MGo] _1 > 

= G + < GoMGqMGo > 

+ < GoMGqMGoMGqMGo > + . . . (18) 

The Feynman graph expansion is shown in fig.(l). One can easily convince oneself that the 
above formulation of the problem amounts to summing all spin-fluctuation exchange graphs. 
Note that no fermion loops are included. The second order dressing of the spin-fluctuation 
propagator and higher order RPA-like diagrams can be absorbed in the definition of the 
magnetic response function x(q). Higher order fermion loops, such as the fourth order term 
with four spin-fluctuation lines (non-linear magnetic response), etc.. are not included in the 
above formulation of the problem. These are the terms that give rise to the "fermion sign 
problem" in determinant Monte Carlo simulations. Some of these diagrams are shown in 
fig. (2). The quartic term should in principle be introduced in the probability distribution of 
the vector field M(q) when the coefficient x _1 (q) °f the quadratic term becomes very small, 
i.e. near the phase transition. In this paper, however, we shall only consider a Gaussian 
distribution of the fields M(q). 
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The above formulation immediately suggests a non-perturbative approach to the prob- 
lem: instead of expanding the matrix inverse and doing the integrals analytically, one can 
sample the functional integral by the Monte Carlo method and numerically invert the Green's 
function matrix in Eq. (|T5|) . The columns of the inverse are calculated by solving the appro- 
priate linear system of equations by means of a preconditioned conjugate gradient algorithm. 
Since the matrix is not hermit ian positive definite, there is no guaranteed convergence for 
the conjugate gradient method. In order to make the algorithm more robust, especially as 
the coupling gets strong, we combine the conjugate gradient with a continuation scheme: 
we invert the matrix 1 — AMGo and gradually switch on the parameter A from zero to one. 
The diagonal of the matrix is used as a preconditioner. 

Since upon doing the averages over all of the field configurations one restores translation 
invariance and other symmetries, the number of linear systems of equations that one must 
solve is much lower than the order of the matrix. 

We have checked that the random fields M were sampled properly by calculating the 
two-point correlation function and comparing to the expected answer Eq. flT2|). A further 
check was to calculate the second order term in the expansion Eq. ( |18D by Monte Carlo 
sampling and compare with the computation of the appropriate Feynman diagram. The 
results agreed within the (small) statistical uncertainty. The conjugate gradient matrix 
inversion was tested on a small system (4 by 4 lattice) at high temperature where it could 
be compared to the result of Gaussian elimination with partial pivoting. 

III. RESULTS 

We shall present results for a doping n=0.8 at a temperature T = 200°K for a couple of 
coupling constants g 2 and compare the results of the computer simulations to lowest order 
perturbation theory (i.e. one spin-fluctuation exchange Feynman graph with bare Green's 
function) and self-consistent one-loop calculations (i.e single spin-fluctuation exchange Feyn- 
man diagram with dressed Green's function). The computations were carried out on a 16 x 16 
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lattice with 40 fermion and 21 boson Matsubara frequencies. With such parameters, the 
Green's function, Eq. (p~5| ) is obtained by inverting the 20480 by 20480 matrix 1 — MG for 
each configuration of the random vector fields M(q). 

The Eliahsberg renormalization factor Z(p, i7rT) for a coupling constant g 2 = 0.1 eV 2 is 
shown in fig. (3). The coupling is intermediate, since there is a difference between the simple 
first-order perturbation theoretic calculation and the self-consistent one-loop result. Note 
that if one tries to calculate the Green's function, Eq. ( |15|) for each configuration of the 
random field M by expanding the denominator in powers of M (using the geometric series), 
one finds that the expansion diverges for most field configurations. We have independently 
checked that the spectral radius of the matrix MG is larger than one, but none of its 
eigenvalues are equal to one, so the inverse in Eq. (|T5l ) is well defined. In this coupling 
constant regime, the self-consistent calculation agrees quite well with the results of the Monte 
Carlo calculation. As expected, it is near the hot spots that the first-order perturbation 
theory and self-consistent calculations deviate from the Monte Carlo results. Away from the 
hot spots, the effective coupling is weak (since there is no appreciable difference between the 
first-order and self-consistent calculations) and perturbation theory works satifactorily. Near 
the hot spots, the Monte Carlo simulation yields a renormalization factor that is slightly 
reduced relative to the self-consistent Eliashberg result. This must be attributed to third 
and higher order processes, since the lowest order vertex correction enhances the effect of 
the interaction for quasiparticles near the hot spots (Chubukov et. al (1997), Monthoux 
(1997)). 

The Eliahsberg renormalization factor Z(p, i7rT) for a coupling constant g 2 = 0.4 eV 2 
is shown in fig. (4). While near cold spots the results of both the first-order perturbation 
theoretic and the self-consistent calculations agree with the Monte Carlo data and each other, 
one notes that near the hot spots the perturbation theoretic calculation does a much better 
job than the Eliashberg self-consistent one. The vertex corrections must thus "undress" the 
Eliashberg Green's function for this value of the coupling constant. Since the results of the 
simulations near the hot spots yield a Z some 20% larger than the Eliashberg calculation, one 
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can conclude that the higher order vertex corrections enhance the spin-fluctuation interaction 
around the antiferromagnetic wavevector when g 2 = 0.4 eV 2 . We have just begun the 
exploration of the parameter space and it will be interesting to see what happens as all of 
the parameters (coupling constant, antiferromagnetic correlation length and characteristic 
spin-fluctuation frequency) are varied. 

IV. OUTLOOK 

We have presented a non-perturbative method for studying the problem of quasiparticles 
coupled to spin-fluctuations and compared our results to first-order perturbation theory and 
one-loop self-consistent calculations. We have shown some of the limitations of both these 
approaches. We shall study the behaviour of the hot quasiparticles as the antiferromagnetic 
length increases and shall look for the pseudo-gap phenomenon seen in ARPES experiments 
(Loeser et al.,1996), (Ding et al.,1996). As one approaches the antiferromagnetic transition, 
it may not be appropriate to consider only Gaussian terms in the action for the spin- 
fluctuations. The method allows the introduction of quartic terms in the action of the 
random magnetic field. The study of the effect of such terms on the normal state properties 
of the system would be of great interest. The present formalism can in principle be extended 
to study other Green's functions. We are presently working on such extensions. 
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FIGURES 
FIG. 1. Spin- fluctuation exchange Feynman graph expansion for the quasiparticle Green's func- 
tion (heavy line) obtained by expanding Eq. (|15|) in powers of the external field M. 

FIG. 2. Some of the diagrams not included in the expansion. The first two RPA-like graphs 
amount to a renormalization of the spin-fluctuation propagator and are already accounted for by 
our use of the experimental %(q, u). The last two diagrams amount to including (M) 4 and (M) 6 
terms in action for the external field M and are not included in this work. 

FIG. 3. Eliashberg renormalization factor Z(p, iu n ) = 1 — ^-ImH(p, icj n ) at T = 200°K as a 
function of momentum p for the lowest Matsubara frequency u n — ttT. The coupling constant is 
g 2 = o.l eV 2 . 

FIG. 4. Eliashberg renormalization factor Z(p, kj n ) = 1 — ^-ImZl(p, icj n ) at T = 200°K as a 
function of momentum p for the lowest Matsubara frequency uj u = 7rT. The coupling constant is 
g 2 = Q.4eV 2 . 
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